--- title: Title keywords: fastai sidebar: home_sidebar nb_path: "nbs/01_HitTest.ipynb" ---
{% raw %} 
51.6º, 179.2º, 38068.6km
Making fig 0.
2017-01-04 20:00:00	data/VAULT_Data/TLE_daily/2017/01/04.tab.gz
Making fig 1.
2017-01-06 03:00:00	data/VAULT_Data/TLE_daily/2017/01/06.tab.gz
Making fig 2.
2017-01-08 18:00:00	data/VAULT_Data/TLE_daily/2017/01/08.tab.gz
Making fig 3.
2017-01-12 05:00:00	data/VAULT_Data/TLE_daily/2017/01/12.tab.gz
Making fig 4.
2017-01-19 09:00:00	data/VAULT_Data/TLE_daily/2017/01/19.tab.gz
Making fig 5.
2017-01-26 18:00:00	data/VAULT_Data/TLE_daily/2017/01/26.tab.gz
Making fig 6.
2017-01-29 18:00:00	data/VAULT_Data/TLE_daily/2017/01/29.tab.gz
{% endraw %}

HitTest

Test a more sensible approach to satellite visibility by using standard astronomy libraries to ask whether a satellite is in view of a ship at a given time.

{% raw %}
{% endraw %} {% raw %}
{% endraw %}

Load the day's TLE file

Create a function to load a single day's TLE file and return parsed datatypes.

{% raw %}
{% endraw %} {% raw %}

load_day_file[source]

load_day_file(_day:datetime, folder:str='data/VAULT_Data/TLE_daily', date_cols=['day_dt', 'tle_dt'], verbose=True)

Look for and load TLE datafile for {_day}.

{% endraw %}

Then test it on a single day.

{% raw %} 
df = load_day_file(datetime(2016, 6, 30))
df.count()
df.head()

2016-06-30 00:00:00	../data/VAULT_Data/TLE_daily/2016/06/30.tab.gz
satellite day_dt day tle_dt tle_ts line1 line2
0 1000 2016-06-30 6026 2016-06-27 11:15:21 1467040521 1 01000U 65008B 16179.46899882 .00000021 0... 2 01000 32.1467 333.7511 0009366 165.3909 194...
1 1000 2016-06-30 6026 2016-06-27 11:15:21 1467040521 1 01000U 65008B 16179.46899882 .00000021 0... 2 01000 32.1467 333.7511 0009366 165.3909 194...
2 1000 2016-06-30 6026 2016-06-27 11:15:21 1467040521 1 01000U 65008B 16179.46899882 .00000021 0... 2 01000 32.1467 333.7511 0009366 165.3909 194...
3 10000 2016-06-30 6026 2016-06-30 10:49:53 1467298193 1 10000U 77034A 16182.45131225 -.00000171 0... 2 10000 15.5820 331.7785 0019081 259.0540 28...
4 10002 2016-06-30 6026 2016-06-28 23:10:32 1467169832 1 10002U 77034C 16180.96565494 -.00000126 0... 2 10002 16.1681 333.0471 0296361 5.9346 0...
{% endraw %} {% raw %} 
df.satellite.value_counts()

29201    106
39694    101
33472     94
28584     87
29203     81
        ... 
333        1
18764      1
29003      1
34004      1
16384      1
Name: satellite, Length: 12152, dtype: int64
{% endraw %}

Drop Dupes

Wait, each satellite should only need one TLE entry. These multiples are all duplicates. Watch.

{% raw %} 
df = df.drop_duplicates()
df.shape

(12152, 7)
{% endraw %}

Memory check

Inspect the resulting dtypes and memory usage.

The parsing was successful. As expected, line1 and line2 are large. Using category doesn't save much because so many rows are unique. The datetime categories are surprisingly large.

{% raw %} 
pd.DataFrame([df.dtypes, df.memory_usage(index=False, deep=True)], index=['Dtype', 'Mem']).T
Dtype Mem
satellite uint16 24304
day_dt datetime64[ns] 97216
day uint16 24304
tle_dt datetime64[ns] 97216
tle_ts uint32 48608
line1 string 1531152
line2 string 1531152
{% endraw %} {% raw %} 
df.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 12152 entries, 0 to 15178
Data columns (total 7 columns):
 #   Column     Non-Null Count  Dtype         
---  ------     --------------  -----         
 0   satellite  12152 non-null  uint16        
 1   day_dt     12152 non-null  datetime64[ns]
 2   day        12152 non-null  uint16        
 3   tle_dt     12152 non-null  datetime64[ns]
 4   tle_ts     12152 non-null  uint32        
 5   line1      12152 non-null  string        
 6   line2      12152 non-null  string        
dtypes: datetime64[ns](2), string(2), uint16(2), uint32(1)
memory usage: 569.6 KB
{% endraw %} {% raw %} 
print(df.iloc[3]["line1"])
print(df.iloc[3]["line2"])

1 01001U 65008A   16178.87975142  .00000008  00000-0  00000+0 0  2425
2 01001  32.1427 144.4001 0016649 349.5981  10.4171  9.90724672860590
{% endraw %}

Skyfield

First, test the basic operation works as expected.

{% raw %} 
51.6º, 179.2º, 38068.6km
{% endraw %} {% raw %}

test_skyfield[source]

test_skyfield()

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
assert (datetime(1971, 6, 1) - datetime(1970, 6, 1)).days == 365
{% endraw %} {% raw %} 
pd.DataFrame([51.5, 189, 2.3], index=['alt','az','days'])
0
alt 51.5
az 189.0
days 2.3
{% endraw %}

Can they see me?

Given at Lat/Lon/Time, what satellites can see me?

We pre-partition the TLE data by Year/Month/Day, so we can quickly load only today's TLE data, and check whether Lat/Lon can see it.

First step: get Alt/Az/dt for each row.

Here we apply the Skyfield.EarthSatellite function to all TLE rows in the dataframe for today.

Benchmark: this takes about 6s on a laptop. @TODO: speed this up by 10x.

Minor glitch: cannot use faster raw=True

In theory apply(..., raw=True) should be faster than default apply. However, it's not working due to:

AssertionError:Number of manager items must equal union of block items> # manager items:7, # tot_items:

  • Possible solution: https://www.nuomiphp.com/eplan/en/254300.html
  • On the other hand, it's an open pandas ticket: https://github.com/pandas-dev/pandas/issues/34822

So I try the default way first. But the except won't work until we track down the block manager issue.

{% raw %}
{% endraw %} {% raw %}

satellite_alt_az_days[source]

satellite_alt_az_days(_t0:datetime, lat:float, lon:float)

Load tracks for day {_t0} and return altitiude, azimuth, and 𝚫t [days] for each row.

{% endraw %}

Execute for a given day

2016-06-30 for starters

{% raw %} 
df_alt_az_days = satellite_alt_az_days(datetime(2016, 6, 30), 45.0, -176.0)

2016-06-30 00:00:00	../data/VAULT_Data/TLE_daily/2016/06/30.tab.gz
{% endraw %} {% raw %} 
df_alt_az_days.count()

altitude    12152
azimuth     12152
days        12152
dtype: int64
{% endraw %} {% raw %} 
df_alt_az_days.head()
altitude azimuth days
0 -42.413819 297.656662 2 days 12:44:39
3 51.646300 179.205373 0 days 10:49:53
4 -32.583116 299.096010 1 days 00:49:28
6 -10.852721 159.700750 3 days 02:53:10
10 -28.903495 101.253621 0 days 04:49:32
{% endraw %}

Second step: Calculate the hit quality

First approximation:

  • It's a hit if the alt > 0 (above the horizon).
  • Smaller time difference -> better quality.

TODO: Kevin, did I capture that logic correctly? I'm confused how a 2-day lag can be "excellent". These aren't days are they?

{% raw %}
{% endraw %} {% raw %}

hit_quality[source]

hit_quality(df_alt_az_days)

Return hit/miss and quality as time proximity.

Parameters

df_alt_az_days: Dataframe returned by satellite_alt_az_days.

Returns

Dataframe with columns ["hit", "miss"]. Each row will have exactly one filled, with a string denoting how recent the pass was, e.g. "excellent", "good", "poor", "stale".

{% endraw %} {% raw %} 
df_hit_quality = hit_quality(df_alt_az_days)
{% endraw %} {% raw %} 
df_hit_quality["hit"].value_counts()

excellent    1490
good           30
poor            1
Name: hit, dtype: int64
{% endraw %} {% raw %} 
df_hit_quality["miss"].value_counts()

excellent    10439
good           178
poor             9
stale            5
Name: miss, dtype: int64
{% endraw %} {% raw %} 
pd.concat([df_hit_quality["hit"].value_counts(), df_hit_quality["miss"].value_counts()], axis=1, sort=False)
hit miss
excellent 1490.0 10439
good 30.0 178
poor 1.0 9
stale NaN 5
{% endraw %} {% raw %} 
df_alt_az_days_visible = df_alt_az_days[df_alt_az_days["altitude"]>0].copy()
{% endraw %} {% raw %} 
df_alt_az_days_visible.count()

altitude    1521
azimuth     1521
days        1521
dtype: int64
{% endraw %} {% raw %} 
df_alt_az_days_visible.head(5)
altitude azimuth days
3 51.646300 179.205373 0 days 10:49:53
17 48.571486 219.021839 0 days 17:13:44
21 35.295619 324.390831 0 days 12:34:30
22 67.176418 358.140393 1 days 15:57:40
25 2.153692 266.529239 0 days 09:03:34
{% endraw %}

Visualize the results

Generate a polar alt/az plot of the qualifying satellites

  • Excellent = blue
  • Good = red
  • Else = yellow

TODO: Is "0" here 0 altitude? That would be on the horizon, which is counter-intuitive.

Note the band of satellites at southern bearings -- this ship was in the Northern hemisphere.

{% raw %}
{% endraw %} {% raw %}

viz[source]

viz(df, show=True, size0=1)

Polar plots a df_alt_az_days_visible dataframe. Dataframe must have: color, days, altitude, azimuth. Returns a Plotly Express polar plot figure. If show=True, also displays it here. size0 is the smallest marker size (used for best hits)

{% endraw %} {% raw %} 
fig = viz(df_alt_az_days_visible)
{% endraw %}

David says try these from 2017:

  • 4-jan-2017, 8pm
  • 6 jan 3am
  • 8 jan 6pm
  • 12 jan 5am
  • 19 jan 9am
  • 26 jan 6pm
  • 29 jan 6pm
{% raw %} 
Making fig 0.
2017-01-04 20:00:00	../data/VAULT_Data/TLE_daily/2017/01/04.tab.gz
Making fig 1.
2017-01-06 03:00:00	../data/VAULT_Data/TLE_daily/2017/01/06.tab.gz
Making fig 2.
2017-01-08 18:00:00	../data/VAULT_Data/TLE_daily/2017/01/08.tab.gz
Making fig 3.
2017-01-12 05:00:00	../data/VAULT_Data/TLE_daily/2017/01/12.tab.gz
Making fig 4.
2017-01-19 09:00:00	../data/VAULT_Data/TLE_daily/2017/01/19.tab.gz
Making fig 5.
2017-01-26 18:00:00	../data/VAULT_Data/TLE_daily/2017/01/26.tab.gz
Making fig 6.
2017-01-29 18:00:00	../data/VAULT_Data/TLE_daily/2017/01/29.tab.gz
{% endraw %} {% raw %}
{% endraw %}